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ABSTRACT 

We consider the effect of a supernova (SN) expiosion in a very massive binary that is 
expected to form in a portion of Population III stars with the mass higher than lOOM©. 
In a Population III binary system, a more massive star can result in the formation of a 
BH and a surrounding accretion disc. Such BH accretion could be a significant source 
of the cosmic reionization in the early universe. However, a less massive companion 
star evolves belatedly and eventually undergoes a SN explosion, so that the accretion 
disc around a BH might be blown off in a lifetime of companion star. In this paper, 
we explore the dynamical impact of a SN explosion on an accretion disc around a 
massive BH, and elucidate whether the BH accretion disc is totally demolished or not. 
For the purpose, we perform three-dimensional hydrodynamic simulations of a very 
massive binary system, where we assume a BH of 10'^ A/© that results from a direct 
collapse of a very massive star and a companion star of lOOM© that undergoes a SN 
explosion. We calculate the remaining mass of a BH accretion disc as a function of 
time. As a result, it is found that a significant portion of gas disc can survive through 
three-dimensional geometrical effects even after the SN explosion of a companion star. 
Even if the SN explosion energy is higher by two orders of magnitude than the binding 
energy of gas disc, about a half of disc can be left over. The results imply that the 
Population III BH accretion disc can be a long-lived lunrinous source, and therefore 
could be an important ionizing source in the early universe. 
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1 INTRODUCTION 

The studies on the formation of Population III (Pop III) 
stars have shown that first stars i n the universe are likely to 
be as massive as 100 — lOOOMp IJBromm. Coppi fc Larson! 
1 19991 : I Abel. Bryan, fc NormanI T2OOI 



20021 ) or to form in 



a bimodal initial mass f unction (IMF) with peaks of sev- 
eral IOOM0 and ~ IM0 (JNakamura fc Umemurall200ll ). Re- 
cently, it has been revealed that a significant fraction of 
Pop III stars can be expected to form in b inary systems 
IjSaigoet al.ll2004l ). iHeger fc WooslevI l|2002l ) have studied 
the nucleosynthetic evolution of Pop III stars, and have 
shown that Pop III stars in the mass range of m > 260Mq 
result in direct black hole formation, while the mass range 
of 25M0 <m< I4OM0 leads to the supernova (SN) explo- 
sions. Hence, if there is mass difference in a Pop III binary. 
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there can form a binary system composed of a BH with 
m > 26OM0 and a massive star with 25M0 < m < 14OM0 . 
Also, recent general relativistic simulations on a supermas- 
sive star (m > 10'^ M0) have shown that if a star pos- 
sesses large angular momentum, a massive gas disc is left 
over around a forming black hole llShibata fc Shapiro I2OO2I : 



IShibatal200J : IShapirol200i : ISekiguchi fc Shibatall2004[ ). The 
mass fraction of gas disc is 10~* — 10~^ of the initial stellar 
mass, depending on the equation of state and the degree of 
rotation. This disc can fuel the BH in the later evolutionary 
stage. Even if a disc left over is quite small, the Roche lobe 
overfiow from a com panion star niay al so fuel the BH in a 
close binary system (jLawlor et al.ll2008l ). 

On the other hand, it is suggested that Pop HI BH 
accretion may be an important cl ue for the reio nization of 
the universe (Madau et al. 2004; Ricotti fc Ostrik er 2004t|). 
For a BH of ~ IQT'Mq, the accretion disc can emit a black- 
body radiation with an effective temperature of l(f~^ K, if 
the mass a c cretio n rate is close to the Eddigton rate (e.g. 
iKato et al.l l| 19981 )). When the disc mass is several IOM0, 
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the accretion timescale is several 10^ years. Therefore, the 
accreting BH can be an ultraviolet radiation source in a 
longer timescale than the lifetime of a massive companion 
star. Then, the accreting BH can be a significant cont ributor 
for the early reionization (|Hirose fc Umemurall2007h . How- 
ever, if the accretion disc is blown off by the SN explosion of 
a companion star, the accretion timescale may be too short 
to play a significant role for the reionization. 

In this paper, we explore the disruption and stripping 
of a gas disc around a BH by a SN explosion of a mas- 
sive companion star. For the purpose, we perform three- 
dimensional hydrodynamic simulations with a second-order 
scheme, AUSM-DV. In particular, to resolve the propaga- 
tion of a blast wave and the deformation of a gas disc with 
high accuracy, three-dimensional generalized curvilinear co- 
ordinates are used. Such a three-dimensional simulation on 
the collision of SN blast wave with a BH accretion disc has 
been never performed so far. With the simulations, we study 
the stripping efficiency as a function of the kinetic energy 
of the SN blast wave, by changing the separation of binary 
and the mass of gas disc. Then, the dependence on the ratio 
of the blast wave kinetic energy to the binding energy of gas 
disc is analysed in detail. 

The paper is organized as follows. In Section 2, the de- 
tails of present numerical model and method are presented. 
In Section 3, we show the numerical results for the defor- 
mation and stripping of gas disc around a BH by the col- 
lision with a SN blast wave. Also, the dependence on the 
model parameters is studied, and the resultant remaining 
mass fraction of disc is argued. Section 4 is devoted to con- 
clusions. 



2 NUMERICAL MODEL AND METHOD 

2.1 Model 

We consider a binary system composed of a black hole of 
IOOOM0 and a Pop HI companion star of lOOM©. In the 
present analyses, the model parameters are the separation 
between the BH and the companion star, the mass of gas 
disc surrounding the BH, and the den s ity pr ofile of gas disc. 
Based on the analysis bv lSaigo et al.l (|2004l ). we set the bi- 
nary separation, a, to be 700 AU as a fiducial case, and a half 
or one-fifth of the fiducial separation is also examined. As for 
the mass of gas disc, considering the disc formation around a 
black hole (IShibata fc Shapirdl2002l : IShibata|[20o3 : IShapird 
l2004l : ISekiguchi fc Shibatall2004l ) and also the mass loss from 
a Po p HI star before a SN explosion (Hcgcr & Woo slevI 
|2002| ). we set the disc mass to be SOM©, lOM© or ZMq. 
Regarding the density distribution of gas disc, we assume a 
power-law type distribution like p{r) oc r" on the equato- 
rial plane. The power-law index n is assumed to be n = 0, 
— 1, or —2. The temperature of gas disc is presumed to be 
lO^K, since the gas disc is irradiated by ultraviolet radiation 
from a companion star before the explosion. The rotation of 
the disc is assumed to be Keplerian, because the self-gravity 
is negligible in all the cases considered. In vertical direc- 
tions of the disc, we determine the density distributions by 
solving a hydrostatic balance. Also, taking the tidal effect 
into acc ount, the g a s disc is t runcated at t he tidal radius 
given bv lPaczvnskl (|l977l ) and iBoflinl (|200ll ). The tidal ra- 
dius is ~ 0.7 Roche lobe radius of the BH. For instance, 



the outer radius of the gas disc is set to be 280AU for a 
fiducial case. The space out of the gas disc is filled with 
tenuous and high temperature gas, since we cannot treat 
a vacuum in the present hydrodynamic scheme. In order to 
minimize its dynamical effect, we set the density to be negli- 
gibly small. Here, the density and temperature are assumed 
to be 2.5 X 10~^®g cm~^ and lO^K, respectively. The param- 
eters adopted are summarized in Table [l] 

To simulate a SN e xplosion, we employ a m ethod 



similar to that taken by Kitavama fc Yoshidal (|2005l ). and 
iBromm. Yoshida fc Hernguistl ( 20031 ) . where the explosion 
energy is inserted as the thermal energy of the gas ejected 
from a SN. The ejected gas has no initial velocity and is 
accelerated by the gas pressure. The temperature is as- 
sumed to be 1.0 X 10*K for all cases, and the density and 
extent of hot g as is ad justed so that the explosion en- 
ergy is lO^^erg (jHeger fc Wooslcv 200^1). If the separation 
is shorter, the density of hot gas is higher and the extent is 
smaller. For model Al, the density of the hot gas is set to 
be psN = 1.0 X 10~^^g cm""^ and the extent is to be 50AU. 
We tested the cases of different density and extent and con- 
firmed that the impact of a SN explosion on the gas disc 
does not have a strong dependence on the choices, but is 
basically determined by the explosion energy. Also, we reg- 
ulate the duration of SN explosio n so that the total ma ss 
ejected by SN explosion is SOM© (JHeger fc WooslevI [2003 ). 

2.2 Numerical Method 

We solve the three-dimensional Euler equation for inviscid 
gas. Basic equations can be written in the conservative form 
as 



dQ_ dE_ 
dt dx 



dP 
dy 



dG 



+ H = 0, 



(1) 



where Q, E, F, G and H are respectively the following vec- 
tors, 
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(e.g. IVinokuilll974l ). where p is the density of gas, u, v and 
w are respectively the x-, y- and z- components of velocities, 
and fx, fy and /z are the x-, y- and z- components of gravity 
force by a black hole. As the gravity force, we take into 
account only the effect of black hole, since it is a major 
component. The pressure p and the total energy per unit 
volume e are related by the equation of state, 

P^{^~l)[e^^{u^+v^+w^)Y (3) 

where 7 is the ratio of specific heats of gas, which is 7 = 5/3 
here. 

The Reynolds number based on the molecular viscosity 
is very large for the present case. Therefore, we have ne- 



(2) 
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glected the effect of viscosity (and only considered the effect 
of ram pressure), i.e. solved the Euler equations. However, 
we should note that possible turbulence would provide a 
tendency to reduce the remaining disc mass or to extinguish 
the disc. Solving such case acco unting for turbulen t eddy 
viscosity will be our future work (JTakeda et al.lll985h . 

To solve the governing equat ion ([T]) numerically, AUSM- 
DV scheme (jWada fc Lioru|l994[ ) is employed. To properly 
adjust the boundary and solve the flow with high accuracy, 
we use the three-dimensional generalized curvilinear coordi- 
nates as shown in Fig. [l] Because we assume the system is 
in a plane-symmetric with respect to the equatorial plane, 
we perform simulations in the computational domain that 
is composed of a hemisphere and an adjacent protuberant 
region. The BH and the SN explosion are placed on the equa- 
torial plane, where the BH is located at the center of hemi- 
sphere, and the SN explosion occurs at the center of pro- 
tuberant region. In Fig. [l] only a half of the computational 
domain is shown after the domain is cut by a plane perpen- 
dicular to the equatorial plane that connects the BH and the 
SN explosion. The number of grid points is 51 x 101 x 31. 
We treat a black hole as a small central hole with absorbing 
boundary, whose radius is 20 AU for the case of a = 700 AU. 
The gas disc is set up around the small hole, and the in- 
ner radius of the disc is assumed to be 50 AU. The grids are 
generated with higher resolution near the black hole and the 
equatorial plane of the gas disc. Such coordinates allow us to 
treat accurately the flow of the gas disc. The computations 
are carried out only in the upper half domain relative to 
the equatorial plane, because we assume that the system is 
symmetric about this plane. The leftward meshes are built 
up to represent a SN explosion. A free boundary condition 
on outer and inner boundaries is adopted in these meshes. 



3 RESULTS 

When a SN explosion produces a supersonic flow with the 
mass of 5OM0 and the velocity of 4000 km/sec, the kinetic 
energy of the blast wave is ~ 1.1 x 10*" erg for the solid 
angle of a disc viewed from the SN explosion. In the fidu- 
cial case, the binding energy of gas disc is evaluated to be 
~ 3.1 X 10** erg. Hence, the kinetic energy of a SN blast 
wave is much larger than the disc binding energy. Then, the 
collision of the blast wave with a disc is anticipated to be 
devastating. However, the disc is deformed by the collision, 
and also streamlines are bent significantly. Therefore, it is 
unclear whether the whole of disc gas is stripped out by the 
collision of the SN blast wave. 



3.1 Disc Deformation and Stripping 

Figure [2] shows the numerical results at early stages for 
model Al. The density distributions and velocity fields are 
shown in the x — y plane with z = Q {lower panels) and in 
x — z plane with y = Q [upper panels). Figure [2{ a) shows the 
initial state, where high-density hot gas is ejected at 700AU 
from a BH. The hot gas expands almost radially (Fig.[2jb)) 
and rushes toward the gas disc (Fig. [5Jc)). At ^ 0.47 years, 
as seen in Fig. [IJd), the bulk of fiow is bent to go through 
the disc. In this stage, some of gas disc is stripped out, but 



much of the blast wave energy can escape owing to the bend- 
ing of gas flow. On the other hand, the fiow directed to the 
disc center is strongly decelerated by the collision with the 
disc, and begins to deform the disc (Fig. [2je)(f)). 

Figure[3]shows the deformation and stripping of the disc 
at later stages. The figure zooms in the disc regions. The 
rotation period at the outer edge of gas disc is ~ 148 years. 
Fig. [3] (a) is the initial state again. As seen in Figures [3] (b) 
and (c), the blast wave pushes the left-side edge of gas disc 
and produces a dent within a few years, because the time- 
scale of the SN explosion is much shorter than the rotation 
period of the outer edge of gas disc. The rotation of disc 
transfers this dent, resulting in the overall deformation of 
the disc (Fig.|3](d)). Through this deformation, the rotation 
balance breaks down, and resultantly some of gas is blow 
out by centrifugal force, as seen in Figure |3] (e) . Eventually, 
the gas disc is settled in a quasi-steady state within about 
100 years (Fig.[3](f)). As a result, roughly 70% of an original 
disc around the BH is left over. In Figure |31 the deformation 
and stripping processes are shown with three-dimensional 
volume rendering visualization. The survival of the gas disc 
after the heaving and ruffling by the blast wave collision is 
clearly shown. 

Previously, I Wheeler et al.l l|l975h argued that in a bi- 
nary system the mass stripping and ablation from a star 
by the impinging balastwave is well expressed by a non- 
dimensional parameter ^, which is deflned as 



^^ 1 MsN R^ /WSN _ ^ 

~ 4 Mc a2 V t;cs 



(4) 



where Msn is the mass of gas expelled by a SN, Mc and 
R are respectively the mass and radius of object affected 
by blast wave, ^sn is the typical speed of the blast wave, 
VcB is the escape velocity. They obtained the mass fraction 
ejected by the stripping and ablation as a function of ^ . We 
can evaluate ^ in the present system, where AfsN = 5OM0, 
Mc = Mdisc = 30Mq, R = 280AU, a = 700AU and usn = 
4000 km/sec. v^a is approximated to be 



/2GMBHy/' 
{ R J ■' 



(5) 



where Mbh = lOOOM© . Then, we flnd * ~ 3.28 in our model 
Al. According to the criterion for a polytropic star with 
n = 3 given bv IWheeler et al.l ((19751), it is predicted that 
about a half of mass is ejected in the case of ^ ~ 3. In the 
present simulation, about 30% of mass is ejected eventually. 
Therefore, we can conclude that ^ is a fairly good measure 
for the mass ejection even in a disk system. However, in the 
disk system, a little more mass remains than the prediction 
for a spherical star. In the collision of blast wave with a 
gas disk, the streamline is more strongly bent above and 
below the disc. Hence, the momentum transferred to the 
gas in a disk is reduced, and more mass can be left over. 
Therefore, the difference from the spherical prediction can 
be reasonably understood in terms of the geometrical effect. 



3.2 Dependence on Binary Separation 

In order to investigate the dependence on the binary sep- 
aration, we perform additional two simulations with a = 
350 AU (model Bl), which is a half separation of model Al, 
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and with a = 140 AU (model B2), which is a one- fifth sep- 
aration of model Al. We assume the mass of the disc to be 
the same as model Al. Therefore, the disc radius becomes 
smaller and the density is higher for these models. Also, the 
energy of the SN explosion, the temperature, and the total 
mass of ejected gas to be the same as model Al. The rota- 
tion period of the outer edge of gas disc for model Bl and B2 
is ~ 52 years and ~ 13 years, respectively. Figure [5] shows 
the density distributions and the velocity fields in the final 
quasi-steady state, which is reached at t = 38.19 years and 
11.95 years for models Bl and B2, respectively. Hence, the 
time to reach the final state is shorter for the shorter separa- 
tion. Interestingly, it is found that more mass is left over for 
shorter binary separation. Roughly 80 % of disc survives for 
model Bl, and nearly 85 % of disc does for model B2. This 
can be understood in terms of the binding energy Eh of the 
disc relative to the kinetic energy Ek of the SN blast wave. 
The ratio Eb/Eu is - 0.1 and ~ 0.2 for models Bl and B2, 
respectively, while Eb/Ek ~ 0.03 for model A. Thus, Eb/Ek 
is thought to be a key physical quantity that determines the 
remaining mass of disc. This point is argued again later. 



3.3 Dependence on Density Distribution 

To check whether the deformation and stripping of disc de- 
pends on the density distribution of gas disc, simulations 
for discs with different density distributions are performed. 
One is a flat density distribution model with n = on the 
equatorial plane (model CI), and the other is a centrally 
concentrated model with n = — 2 (model C2). The mass of 
gas disc is the same as model Al. Figure |6] shows the re- 
sults for a flat density distribution disc (model CI). The 
snapshots of density distributions and velocity fields at t = 
0.00, 10.32, and 110.16 years are presented. The results look 
similar those shown in Figure [3] at the same evolutionary 
stages. Actually, the remaining mass is not changed sigrufi- 
cantly. Figure[7]show the results for a centrally-concentrated 
disc (model C2). The snapshots of density distributions and 
velocity fields at t = 0.00, 10.37, and 105.69 years are pre- 
sented. It is found that, for model C2, the infiuence by the 
blast wave is slightly weaker in the central regions of the 
disc (Fig. [7](b)), compared to model CI (Fig. [6] (b)). But, 
the remaining mass increases by only several %. Therefore, 
we can conclude that the density distributions of disc do not 
infiuence the final mass of disc strongly. 



3.4 Remaining Mass 

In Figure [S] the time sequences of the disc mass are summa- 
rized by the ratio of the remaining disc mass to the disc mass 
without a supernova explosion of companion star. For model 
Al, CI and C2, the disc is settled in quasi-steady states in 60 
years, and ~ 70% of disc gas remains, almost regardless of 
density distributions of disc. For model Bl (a = 350AU), a 
quEisi-steady state is attained in 20 years and '^ 80% of disc 
gas remains. For model B2 (a — 140AU), the disc reaches a 
quasi-steady state in 10 years and ~ 85% of disc gas remains. 
As argued in Section 13.11 the remaining amount of gas 
can be understood in terms of the ratio of Eb to Ek, where 
Eb is the binding energy of the gas disc around the BH 
and Ek is the kinetic energy of the blast wave taking into 



account the solid angle of the disc viewed from the SN explo- 
sion point. Eb is evaluated by integrating potential energy 
of the gas around the BH. Figure [§] shows the remaining 
mass fraction against the Eb/Ek ratio. In this diagram, two 
more simulations are added with changing the disc mass: 
model Dl has the disc mass of IOM0, which is one-third 
of that in model A, and model D2 has the disc mass of 
SMq, which is one-tenth of that in model A. Obviously, the 
remaining mass fraction correlates positively with Eb/Ek- 
Thus, a key physical parameter that determines the remain- 
ing mass fraction is the Eb/Ek ratio. Interestingly, even if 
the binding energy is less than 1% of the kinetic energy of 
SN blast wave, a considerable amount of mass remains. This 
is a three-dimensional effect of deformation and stripping, 
which is demonstrated in Figure]?] 



4 CONCLUSIONS 

We have explored the dynamical impact of SN explosion in 
a very massive binary system. We have simulated the de- 
formation and stripping of a gas disc around a BH, using 
three-dimensional generalized curvilinear coordinates. Also, 
the dependence on the binary separation, the density dis- 
tribution of disc, and the mass of disc have been studied. 
As a result, we have found that a SN blast wave is not so 
devastating that a disc around a BH is totally evaporated. It 
has turned out that the remaining mass fraction is basically 
determined by the ratio of the disc binding energy to the 
SN kinetic energy. As a matter of importance, even if the 
binding energy of gas disc is much smaller than the kinetic 
energy of blast wave, a significant amount of mass can re- 
main. For instance, when the binding energy is about 1% of 
the kinetic energy of SN blast wave, about a half of mass 
is left over. We have found that the three dimensionality of 
the disc deformation is essential to make a disc around a BH 
survive. 

The present results imply that in a massive binary sys- 
tem like a Pop HI binary, the BH accretion activity can 
be long-lived even after a massive companion star explodes. 
Such BH accretion can emit ultraviolet radiation, and there- 
fore a Pop HI BH accretion disc can be an important can- 
didate as a reionization source in the early universe. 
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Table 1. Adopted parameters for model calculations. In this table, a represents the binary separation, M^iac does the mass of gas disc, 
and n does the power-law index of density distributions on the equatorial plane. 



Model a [AU] Mjis, 



Remarks 



Al 


700 


30 


-1 


Bl 


350 


30 


-1 


B2 


140 


30 


-1 


CI 


700 


30 





C2 


700 


30 


-2 


Dl 


700 


10 


-1 


D2 


700 


3 


-1 



Fiducial model 

Half of the separation in model Al 

One-fifth of the separation in model Al 

Flat density distribution disc 

Centrally concentrated density distribution disc 

One-third of the disc mass in model Al 

One-tenth of the disc mass in model Al 




Figure 1. The three-dimensional generalized curvilinear coordinate system used in this paper. The bottom plane is the equatorial plane. 
The center of the gas disc, i.e. the BH, is placed at the center of a rightward hemisphere, and the SN explosion occurs at the center of an 
adjacent protuberant region. A half of grids in the computational domain are shown, after the whole domain is cut by a plane including 
the BH and the SN explosion point perpendicular to the equatorial plane. The number of grid points is 51 X 101 X 31. 
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Figure 2. Early stages of the interaction of a SN blast wave with a disc around a BH for model Al. Snapshots show density distributions 
and velocity fields in the x — y plane with 2 = [lower panels)^ and in a; — 2: plane with y = {upper panels). 
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Figure 3. Later stages of the interaction of a SN blast wave with a disc around a BH for model Al, where the deformation and stripping 
processes of the gas disc are shown. The density distributions and the velocity fields are shown with zooming in disc regions. 
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Figure 4. 3D volume rendering visualization of density distributions with velocity fields for model Al. Figures (a), (b), (c), (d), (e), 
(f), (g), (h) and (i) represent the snapshots at i = 0.00, 0.63, 0.78, 2.90, 11.17, 35.72, 53.98, 72.98, and 110.24 yrs, respectively. 
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Figure 5. The density distributions and the velocity fields in the final state for shorter separation models; model Bl is a model with 
a = 350 AU (left panel) and model B2 is a model with a = 140 AU {right panel). The snapshots at i = 38.19 years and 11.95 years are 
presented for model Bl and B2, respectively. 
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Figure 6. Same as Fig.[3]but for a flat density distribution model with n = (model CI). The snapshots at £ = 0.00, 10.32, and 110.16 
years are presented. 
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Figure 7. Same as Fig. [3] but for a central- concentration density distribution model with n = —2 (model C2). The snapshots at 
t = 0.00, 10.37, and 105.69 years are presented. 
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Figure 8. The time variation of the remaining mass of a disc affected by SN explosion, where the ratio of remaining mass to the mass 
of an undisturbed disc is shown. The horizontal axis is time in units of years. A solid line, a dashed line, a dash-dotted line, a dotted 
line and a log-dashed line show model Al, Bl, B2, CI and C2, respectively. Note that the ratio for each case starts exactly from unity. 
But, the strong inflow by supernova explosion causes the transient rapid growth of ratio by about 10% at the very early stage. 
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Figure 9. The resultant remaining mass fraction is shown against the ratio of E^ to _Efe, where _Ef, is the binding energy of the gas disc 
around a BH and E}^ is the kinetic energy of blast wave taking into account the solid angle of the disc viewed from the SN explosion 
point. The names near the symbols arc model names. 



